home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / newmat03.lha / newmat03 / fft.cxx < prev    next >
C/C++ Source or Header  |  1993-08-08  |  4KB  |  137 lines
  1. //$$ fft.cxx                         Fast fourier transform
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5.  
  6. #define WANT_MATH
  7.  
  8. #include "include.hxx"
  9.  
  10. #include "newmatap.hxx"
  11.  
  12.  
  13. static void cossin(int n, int d, real& c, real& s)
  14. // calculate cos(twopi*n/d) and sin(twopi*n/d)
  15. // minimise roundoff error
  16. {
  17.    long n4 = n * 4; int sector = (int)floor( (real)n4 / (real)d + 0.5 );
  18.    n4 -= sector * d;
  19.    if (sector < 0) sector = 3 - (3 - sector) % 4; else sector %= 4;
  20.    real ratio = 1.5707963267948966192 * (real)n4 / (real)d;
  21.  
  22.    switch (sector)
  23.    {
  24.    case 0: c =  cos(ratio); s =  sin(ratio); break;
  25.    case 1: c = -sin(ratio); s =  cos(ratio); break;
  26.    case 2: c = -cos(ratio); s = -sin(ratio); break;
  27.    case 3: c =  sin(ratio); s = -cos(ratio); break;
  28.    }
  29. }
  30.  
  31. static void fftstep(ColumnVector& A, ColumnVector& B, ColumnVector& X,
  32.    ColumnVector& Y, int after, int now, int before)
  33. {
  34.    // const real twopi = 6.2831853071795864769;
  35.    const int gamma = after * before;  const int delta = now * after;
  36.    // const real angle = twopi / delta;  real temp;
  37.    // real r_omega = cos(angle);  real i_omega = -sin(angle);
  38.    real r_arg = 1.0;  real i_arg = 0.0;
  39.    real* x = X.Store();  real* y = Y.Store();   // pointers to array storage
  40.    const int m = A.Nrows() - gamma;
  41.  
  42.    for (int j = 0; j < now; j++)
  43.    {
  44.       real* a = A.Store(); real* b = B.Store(); // pointers to array storage
  45.       real* x1 = x; real* y1 = y; x += after; y += after;
  46.       for (int ia = 0; ia < after; ia++)
  47.       {
  48.      // generate sins & cosines explicitly rather than iteratively
  49.      // for more accuracy; but slower
  50.      cossin(-(j*after+ia), delta, r_arg, i_arg);
  51.  
  52.      real* a1 = a++; real* b1 = b++; real* x2 = x1++; real* y2 = y1++;
  53.      if (now==2)
  54.      {
  55.         int ib = before; while (ib--)
  56.         {
  57.            real* a2 = m + a1; real* b2 = m + b1; a1 += after; b1 += after;
  58.            real r_value = *a2; real i_value = *b2;
  59.            *x2 = r_value * r_arg - i_value * i_arg + *(a2-gamma);
  60.            *y2 = r_value * i_arg + i_value * r_arg + *(b2-gamma);
  61.            x2 += delta; y2 += delta;
  62.         }
  63.      }
  64.      else
  65.      {
  66.         int ib = before; while (ib--)
  67.         {
  68.            real* a2 = m + a1; real* b2 = m + b1; a1 += after; b1 += after;
  69.            real r_value = *a2; real i_value = *b2;
  70.            int in = now-1; while (in--)
  71.            {
  72.           // it should be possible to make this faster
  73.           // hand code for now = 2,3,4,5,8
  74.           // use symmetry to halve number of operations
  75.           a2 -= gamma; b2 -= gamma;  real temp = r_value;
  76.           r_value = r_value * r_arg - i_value * i_arg + *a2;
  77.           i_value = temp    * i_arg + i_value * r_arg + *b2;
  78.            }
  79.            *x2 = r_value; *y2 = i_value;   x2 += delta; y2 += delta;
  80.         }
  81.      }
  82.  
  83.          // temp = r_arg;
  84.          // r_arg = r_arg * r_omega - i_arg * i_omega;
  85.          // i_arg = temp  * i_omega + i_arg * r_omega;
  86.  
  87.       }
  88.    }
  89. }
  90.  
  91.  
  92. void FFT(const ColumnVector& U, const ColumnVector& V,
  93.    ColumnVector& X, ColumnVector& Y)
  94. {
  95.    // from Carl de Boor (1980), Siam J Sci Stat Comput, 1 173-8
  96.    const int n = U.Nrows();                     // length of arrays
  97.    if (n != V.Nrows()) MatrixError("FFT - vector lengths unequal");
  98.    if (n == 0) MatrixError("FFT - vector length zero");
  99. #ifdef __ZTC__
  100.    ColumnVector A = U.c(); ColumnVector B = V.c();
  101. #else
  102.    ColumnVector A = U; ColumnVector B = V;
  103. #endif
  104.    X.ReDimension(n); Y.ReDimension(n);
  105.    const int nextmx = 8;
  106. #ifndef ATandT
  107.    int prime[8] = { 2,3,5,7,11,13,17,19 };
  108. #else
  109.    int prime[8];
  110.    prime[0]=2; prime[1]=3; prime[2]=5; prime[3]=7;
  111.    prime[4]=11; prime[5]=13; prime[6]=17; prime[7]=19;
  112. #endif
  113.    int after = 1; int before = n; int next = 0; BOOL inzee = TRUE;
  114.  
  115.    do
  116.    {
  117.       int now, b1;
  118.       for (;;)
  119.       {
  120.      if (next < nextmx) now = prime[next];
  121.      b1 = before / now;  if (b1 * now == before) break;
  122.      next++; now += 2;
  123.       }
  124.       before = b1;
  125.  
  126.       if (inzee) fftstep(A, B, X, Y, after, now, before);
  127.       else fftstep(X, Y, A, B, after, now, before);
  128.  
  129.       inzee = !inzee; after *= now;
  130.    }
  131.    while (before != 1);
  132.  
  133.    if (inzee) { A.Release(); X = A; B.Release(); Y = B; }
  134. }
  135.  
  136.  
  137.